!***************************************************************************************************************************
! Program F_Exporter will export data from G-Range, indicated by parameters passed to it in the command line.
! Six parameters are used, they are:
!   1.  The pathway to the G-Range output file storing the data to be exported (e.g., C:\GRange\G_Range\Output)
!   2.  The name of the G-Range output file storing the data to be exported (e.g., pot_evap.gof)
!   3.  The pathway and root portion of a file name for the data that are to be exported (e.g., C:\GRange2\Exp_Out\PET)
!   4.  The first year to be exported (e.g., 2004)
!   5.  The last year to be exported (e.g., 2006)
!   6.  The number of layers in the variable to be exported (e.g., 1, 3, 6).  Some G-Range variables have a single layer 
!           (e.g., bare_ground) and some have multiple layers (e.g., facets, with herbs, shrubs, and trees).
! Output is as GRIDASCII files.
!
! Programmer:  RB Boone   Date:  12/21/2013  Modified Date:  12/28/2013
!***************************************************************************************************************************
program main
  use Parms
  
  echo_on = .FALSE.
  call get_parms
  call load_range_cells
  call load_globe
 
  do year=start_year, end_year
    do month=1, 12 
      call get_month_string   ! Set a variable to a two-character string.  
      do layer=1,layers
        write(*,*) "Exporting layer", layer,"of month",month,"of year",year
        call export_month
      end do      
    end do
  end do

end program main


!***************************************************************************************************************************
! Load the rangeland cells, the cells modeled by G-Range
!***************************************************************************************************************************
subroutine load_range_cells
  use Parms
  integer ::  i
  integer ::  iLower_X, iLower_Y, iUpper_X, iUpper_Y

  open(1, FILE=trim(in_path)//"Rng_Data.gof", ACTION='READ', IOSTAT=ioerr)
   if (ioerr == 0) then
     if (echo_on .eq. .TRUE.) write(*,*) 'Opening the Rng_Data.gof file'
     ! Read the header information
     read(1,*) r_cell_cnt, x_dim, y_dim, lower_x, lower_y, upper_x, upper_y, cell_size, first_year, last_year
     if (echo_on .eq. .TRUE.) then
       write(*,*) 'Range cell count: ', r_cell_cnt
       write(*,*) 'X dimension:      ', x_dim
       write(*,*) 'Y dimension:      ', y_dim
       write(*,*) 'Lower X:          ', lower_x
       write(*,*) 'Lower Y:          ', lower_y
       write(*,*) 'Upper X:          ', upper_x
       write(*,*) 'Upper Y:          ', upper_y
       write(*,*) 'Cell size:        ', cell_size
       write(*,*) 'First year:       ', first_year
       write(*,*) 'Last year:        ', last_year
     end if
     iLower_Y = ((upper_y - lower_y) / 2.0) / cell_size 
     iLower_Y = sign(iLower_Y, -1)
     iUpper_Y = ((upper_y - lower_y) / 2.0) / cell_size
     iLower_X = ((upper_x - lower_x) / 2.0) / cell_size 
     iLower_X = sign(iLower_X, -1)
     iUpper_X = ((upper_x - lower_x) / 2.0) / cell_size
     do i=1, r_cell_cnt
       read(1,*) r_cells(i)%cell_index, r_cells(i)%zone, r_cells(i)%x, r_cells(i)%y, r_cells(i)%range_type
       r_cells(i)%x = r_cells(i)%x + iLower_X
       r_cells(i)%y = r_cells(i)%y + iLower_Y
       if (echo_on .eq. .TRUE.) then
         write(*,*) 'Cells, Zone, X, Y, Range_Type: ', &
         r_cells(i)%cell_index, r_cells(i)%zone, r_cells(i)%x, r_cells(i)%y, r_cells(i)%range_type
       end if     
     end do
   else
     write(*,*) 'The Rng_Data.gof file cannot be opened.'
     stop
   end if
   close(1)

end subroutine


!***************************************************************************************************************************
! Load the globe, used to identify land and ocean cells
!***************************************************************************************************************************
subroutine load_globe
  use Parms
  integer ::  xer, yer
  integer ::  iLower_Y, iUpper_Y, iLower_X, iUpper_X
  character(4350)    ::  globe_line             ! More than is needed for (-2160:2160)

  open(2, FILE=trim(in_path)//"Globe.gof", ACTION='READ', IOSTAT=ioerr)
  if (ioerr == 0) then
    if (echo_on .eq. .TRUE.) write(*,*) 'Opening the Globe.gof file'
    ! Required to change the program to accomodate real cell sizes.  The logic worked for integers (1 degree) but fails otherwise.
    iLower_Y = ((upper_y - lower_y) / 2.0) / cell_size 
    iLower_Y = sign(iLower_Y, -1)
    iUpper_Y = ((upper_y - lower_y) / 2.0) / cell_size
    iLower_X = ((upper_x - lower_x) / 2.0) / cell_size 
    iLower_X = sign(iLower_X, -1)
    iUpper_X = ((upper_x - lower_x) / 2.0) / cell_size
  
    do yer=iLower_Y,iUpper_Y - 1                  ! Adjusting as needed to align map elemenets (was - 0, - 1).  Experimental.
      read(2,*) globe_line
      do xer=iLower_X + 1, iUpper_X                  ! Adjusting as needed to align map elemenets (was - 0, - 1).  Experimental.
        if (trim(globe_line(xer-iLower_X + 1:xer-iLower_X + 1)) .eq. "T") then        
          globe(xer,yer) = .TRUE.
        else
          globe(xer,yer) = .FALSE.
        end if
      end do
    end do
    if (echo_on .eq. .TRUE.) then
      do yer=iUpper_Y -1, iLower_Y, -1
        globe_line = ""
        do xer=iLower_X, iUpper_X - 1
          if (globe(xer - iLower_X,yer) .eq. .TRUE.) then
            globe_line(xer-iLower_X + 1:xer-iLower_X + 1) = "T"
          else
            globe_line(xer-iLower_X + 1:xer-iLower_X + 1) = "F"        
          end if
        end do
        write(*,*) 'ECHOED: ',trim(globe_line)
      end do
    end if
  else
    write(*,*) 'The Globe.gof file cannot be opened.'
    stop
  end if
  close(2)
end subroutine


!***************************************************************************************************************************
! Use unix libraries to get command line arguments.  If anything is missing print a warning and stop
!***************************************************************************************************************************
subroutine get_parms
   use unix_library
   use Parms
 
   character(len=100) :: arg
   logical            :: stopper
   integer            :: i
    
   stopper = .FALSE.

   ! PARM 1
   call getarg(1, arg)                          ! The pathway to the G-Range files
   if ( len(trim(arg)) .eq. 0 .or. ( trim(arg) .eq. "-h" .or. trim(arg) .eq. "help" ) ) then
     write(*,*) 'The first parameter should be the pathway to the G-Range files.'
     stopper = .TRUE.
   else 
     in_path = arg                               ! Not checking valididty of path of path for now.
     i = len(trim(in_path))
     ! *** NOTE: As long as this is compiled for WINDOWS, this will be ok.   If compiled on UNIX,
     ! ***       an alternative separator may be tested for.
     if (in_path(i:i) .ne. "\") then
       in_path = trim(in_path)//"\" 
     end if
   end if
   ! PARM 2
   call getarg(2, arg)                          ! The name of the file to be processed
   if ( len(trim(arg)) .eq. 0 ) then
     write(*,*) 'The second parameter should be the name of the input file to be processed.'
     stopper = .TRUE.
   else 
    file_name = arg                              ! Not checking the existing of the file for now.
   end if
   ! PARM 3
   call getarg(3, arg)                          ! The pathway and root of the files to be created
   if ( len(trim(arg)) .eq. 0 ) then
     write(*,*) 'The third parameter should be the pathway and root of the files to be created.'
     stopper = .TRUE.
   else 
    out_path_root = arg                          ! Not checking the validity of the path for now.
   end if
   ! PARM 4
   call getarg(4, arg)                           ! The first year to be processed.
   if ( len(trim(arg)) .eq. 0 ) then
     write(*,*) 'The fourth parameter should be the first year to be processed.'
     stopper = .TRUE.
   else 
     read(arg,'(I4)') start_year                 ! Using internal conversion
   end if
   ! PARM 5
   call getarg(5, arg)                           ! The last year to be processed.
   if ( len(trim(arg)) .eq. 0 ) then
     write(*,*) 'The fifth parameter should be the last year to be processed.'
     stopper = .TRUE.
   else 
     read(arg,'(I4)') end_year                  ! Using internal conversion
   end if
   ! PARM 6
   call getarg(6, arg)                          ! The layes to be processed.
   if ( len(trim(arg)) .eq. 0 ) then
     write(*,*) 'The sixth parameter should be number of layers to be processed.'
     stopper = .TRUE.
   else 
     read(arg,'(I2)') layers                     ! Using internal conversion
   end if
   if (stopper .eq. .TRUE.) then
     stop
   else
     write(*,*) " "
     write(*,*) "Exporting data from ", trim(in_path)//trim(file_name), " to files named like ", trim(out_path_root), &
          " starting with year", start_year, "and ending with year", end_year, "and including", layers, " layer(s)."
   end if
end subroutine     


!***************************************************************************************************************************
! The subroutine that does the hard work, exporting the layer requested.
!***************************************************************************************************************************
subroutine export_month
  use Parms
 
  integer      :: iLower_X, iLower_Y, iUpper_X, iUpper_Y, icell
  integer*4    :: rec_num                                    ! WHAT IS THE CORRECT RECORD LENGTH?
  integer      :: pre_years, pre_months, pre_images, pre_img_cnt, xer, yer
  real         :: rng_check, yr_check, mn_check              ! WHAT IS THE CORRECT INTEGER LENGTH?
  real         :: v                                          ! WHAT IS THE CORRECT REAL LENGTH?
  character(4) :: year_string
  character(1) :: layer_string
  
  write(year_string,'(I4)') year                  ! Using internal conversion
  write(layer_string,'(I1)') layer
  
  open(21, FILE=trim(out_path_root)//"_"//year_string//"_"//month_string//"_"//layer_string//".asc", ACTION='WRITE', IOSTAT=ioerr)
  if (ioerr == 0) then
    if (echo_on .eq. .TRUE.) write(*,*) 'Opening the file: ', &
      trim(out_path_root)//"_"//year_string//"_"//month_string//"_"//layer_string//".asc"
    call do_headers 
    ! Required to change the program to accomodate real cell sizes.  The logic worked for integers (1 degree) but fails otherwise.
    iLower_Y = ((upper_y - lower_y) / 2.0) / cell_size 
    iLower_Y = sign(iLower_Y, -1)
    iUpper_Y = ((upper_y - lower_y) / 2.0) / cell_size
    iLower_X = ((upper_x - lower_x) / 2.0) / cell_size 
    iLower_X = sign(iLower_X, -1)
    iUpper_X = ((upper_x - lower_x) / 2.0) / cell_size
  
    do yer=iUpper_Y, iLower_Y - 1, -1              ! Must make adjustments to align rangelands with the other global lands
      do xer=iLower_X, iUpper_X
        tglobe(xer, yer) = -255.0
        if (globe(xer, yer) .eq. .FALSE.) Then
          tglobe(xer, yer) = -9999
        end if
      end do
    end do

    close(1)
    open(1, FILE=trim(in_path)//file_name, ACCESS='DIRECT', RECL=4, IOSTAT=ioerr)
    if (ioerr == 0) then
      ! Year times rangecells times the twelve months, plus the current month times rangecells, plus the rangecells checksum.
      ! A rangecells value is written at the 'top' of each map, plus the Year and the Month  
      pre_years = year - first_year
      pre_months = month - 1
      pre_images = layer - 1

      pre_img_cnt = pre_years * 12 * layers * ( r_cell_cnt + 3 )                 ! Years X 12 X Layers X Range cells + the 3 checking values
      pre_img_cnt = pre_img_cnt + ( pre_months * layers * ( r_cell_cnt + 3 ) )   ! Month X Layers X Range cells
      pre_img_cnt = pre_img_cnt + ( pre_images * ( r_cell_cnt + 3 ) )            ! Layers X Range cells
      rec_num = pre_img_cnt + 1                                                  ! Rec_Num starts at 1
      read(1,REC=rec_num) rng_check
      rec_num = rec_num + 1
      read(1,REC=rec_num) yr_check
      rec_num = rec_num + 1
      read(1,REC=rec_num) mn_check
      rec_num = rec_num + 1

      ! Check the various values
      if (rng_check .ne. r_cell_cnt) then
        write(*,*) "There is disagreement between the number of rangeland cells and what is being read in. "
        write(*,*) "Something is out of sync. Expected  was", r_cell_cnt, "and read-in was", rng_check
        stop
      end if
      if (yr_check .ne. year) then
        write(*,*) "There is disagreement between the current year and the year read-in. "
        write(*,*) "Something is out of sync.  Expected is", year, "and read-in was", yr_check
      end if
      if (mn_check .ne. month) then
        write(*,*) "There is disagreement between the current month and the month read-in. "
        write(*,*) "Something is out of sync.  Expected is", month, "and read-in was", mn_check
      end if

      do icell=1,r_cell_cnt
        read(1,REC=rec_num) v
        rec_num = rec_num + 1
        ! The following offsets were made based on comparisons between the land/ocean layer and data layers
        ! 12/28/2013 Making adjustments to the offsets on the next line to align the results with the Visual Basic version.
        tglobe(r_cells(icell)%x - 1, r_cells(icell)%y - 1) = v 
      end do
      close(1)
  
      ! Produce the gridascii file
      do yer=iLower_Y, (iUpper_Y - 1)
        write(21,*) (tglobe(xer, yer), " ",xer=iLower_X, (iUpper_X -1))
      end do
    else
      write(*,*) 'The file to be read from cannot be opened. '
      write(*,*) 'Stopping while attempting to read: ', trim(in_path)//file_name
      stop
    end if    
  else
    write(*,*) 'An output file cannot be created.  Stopping while attempting to create:', &
      trim(out_path_root)//"_"//year_string//"_"//month_string//"_"//layer_string//".asc"
    stop
  end if
  close(21)
  close(1)

end subroutine
  

!***************************************************************************************************************************
! Assign a two character month identifier using brute force
!***************************************************************************************************************************
subroutine get_month_string
  use Parms

  select case ( month )
    case (1)
      month_string = "01" 
    case (2)
      month_string = "02"
    case (3)
      month_string = "03"
    case (4)
      month_string = "04"
    case (5)
      month_string = "05"
    case (6)
      month_string = "06"
    case (7)
      month_string = "07"
    case (8)
      month_string = "08"
    case (9)
      month_string = "09"
    case (10)
      month_string = "10"
    case (11)
      month_string = "11"
    case (12)
      month_string = "12"
  end select

end subroutine


!***************************************************************************************************************************
! Write the appropriate header for the GRIDASCII file
!***************************************************************************************************************************
subroutine do_headers
  use Parms
  
  ! (Select case does not allow floating point case variables)
  if (cell_size .gt. 0.99 .and. cell_size .lt. 1.01) then
    write(21,'(A)') "ncols         360"
    write(21,'(A)') "nrows         180"
    write(21,'(A)') "xllcorner     -180"
    write(21,'(A)') "yllcorner     -90"
    write(21,'(A)') "cellsize      1.0"
    write(21,'(A)') "NODATA_value  -9999"
  end if
  if (cell_size .gt. 0.49 .and. cell_size .lt. 0.51) then   
    write(21,'(A)') "ncols         720"
    write(21,'(A)') "nrows         360"
    write(21,'(A)') "xllcorner     -180"
    write(21,'(A)') "yllcorner     -90"
    write(21,'(A)') "cellsize      0.5"
    write(21,'(A)') "NODATA_value  -9999"
  end if
  if (cell_size .gt. 0.24 .and. cell_size .lt. 0.26) then
    write(21,'(A)') "ncols         1440"
    write(21,'(A)') "nrows         720"
    write(21,'(A)') "xllcorner     -180"
    write(21,'(A)') "yllcorner     -90"
    write(21,'(A)') "cellsize      0.25"
    write(21,'(A)') "NODATA_value  -9999"
  end if
  if (cell_size .gt. 0.16 .and. cell_size .lt. 0.173) then
    write(21,'(A)') "ncols         2160"
    write(21,'(A)') "nrows         1080"
    write(21,'(A)') "xllcorner     -180"
    write(21,'(A)') "yllcorner     -90"
    write(21,'(A)') "cellsize      0.166666666667"
    write(21,'(A)') "NODATA_value  -9999"
  end if
  if (cell_size .gt. 0.09 .and. cell_size .lt. 0.11) then
    write(21,'(A)') "ncols         3600"
    write(21,'(A)') "nrows         1800"
    write(21,'(A)') "xllcorner     -180"
    write(21,'(A)') "yllcorner     -90"
    write(21,'(A)') "cellsize      0.1"
    write(21,'(A)') "NODATA_value  -9999"
  end if
  if (cell_size .gt. 0.08 .and. cell_size .lt. 0.086) then
    write(21,'(A)') "ncols         4320"
    write(21,'(A)') "nrows         2160"
    write(21,'(A)') "xllcorner     -180"
    write(21,'(A)') "yllcorner     -90"
    write(21,'(A)') "cellsize      0.083333333334"
    write(21,'(A)') "NODATA_value  -9999"
  end if
end subroutine
